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0\ ; 1 Introduction 



Long baseline interferometers to detect gravitational waves such as TAMA300 §, GEO600 @, VIRGO 
Pi, LIGO 0] will be in operation by the end of this century or the beginning of the 21st century. One 
of the most important sources of gravitational waves for such detectors are coalescing binary neutron 
\Q ■ stars. PSR1913+16 is the first binary neutron stars found by Hulse and Taylor. The existence of the 

gravitational waves is confirmed by analyzing the arrival time of radio pulses from PSR1913+16. At 
present there are three systems like PSR1913+16 for which coalescing time due to the emission of 
gravitational waves is less than the age of the universe. From these observed systems the event rate 
of the coalescing binary neutron stars is estimated as 10~ 6 ~ 2x 10~ 5 events yr _1 galaxy -1 H ||, 
| From the formation theory of binary neutron stars, the event rate is estimated as 2 x 1CP 5 ~3x 1CP 4 

events yr -1 galaxy -1 [|[ [| |l(J. If the event rate is ~ 2 x 10 -5 , two different methods yield the 
agreement. Since the number density of galaxies is ~ 10 2 /Mpc 3 , we may expect several coalescing 
events/year within 200Mpc (~ 2 x 10 _7 events /Mpc 3 ). The coalescing binary neutron stars is a 
possible central engine of cosmological gamma ray bursts. Coalescing binary neutron stars is also a 
possible site of r-process element production. 

The eccentricity e and the semimajor axis a of binary neutron stars decrease due to the emission 
of gravitational waves, and these quantities are related with each other as 
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j_h ' where eo and do are the initial eccentricity and the semimajor axis, respectively [IllL |12j. If a ~ Rq 

like PSR1913+16, e ~ 10~ 6 for a ~ 500km so that the orbit is almost circular. The merging time 
i mrg due to the emission of the gravitational waves is given by 
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where mj and rri2 are each mass of neutron stars. The frequency of the gravitational waves v GW is 
given by 

: 19Hz (mi+m 2 \^ayi 

V 2.8MP, J V 470km/ v ' 



This is just the frequency range {v GW = 10Hz ~ 1kHz) of long baseline interferometers to detect 
gravitational waves such as TAMA300, GEO600, VIRGO, LIGO. The number of rotation before the 
final merge is given by 
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The merging phase of coalescing binary neutron stars is divided into two stages; l)The last three 
minutes and 2) The last three milliseconds. In the first stage a is much larger than the neutron 
star radius i?NS and the rotation velocity of the binary v r is ~ 0.05c so that each neutron star can be 
considered as a point particle. The post-Newtonian expansion will converge in the first stage. However 
the quite accurate theoretical calculations of the wave form are needed since only the uncertainty of 
1/2635 in the theoretical prediction of the energy loss will cause one rotation ambiguity at last. If the 
accurate theoretical template is obtained |l4|, by making a cross correlation with the observational 
data we may determine each mass, spin, inclination and the distance of binary neutron stars p5| . 

In the second stage a is comparable to i?Ns an d the gravity is so strong that the post-Newtonian 
expansion is not a good approximation and the finite size effect of neutron stars is important. In the 
final collision phase, the general relativistic hydrodynamics are also important. Only 3D numerical 
relativity can study this important final merging phase related to gravitational waves, gamma ray 
burst and production of r-process elements. 



2 Newtonian Simulations 

To clarify the phenomena in the last three milliseconds many numerical simulations have been done 
so far. They are classified into several categories. 

A) Newtonian hydrodynamics without radiation reaction 

In this category, using the Newtonian hydrodynamics, equilibria and stability of the binary neutron 
stars are studied for various initial parameters and equation of states [[l6[ [l7], [fsl E^, ^IL ^2| . 

B) A)+ simple radiation reaction up to contact of neutron stars 

In this category, radiation reaction is included so that the binary spirals in due to the emission of 
the gravitational waves. However after the contact the radiation reaction is switched off p5[ . 

C) A) + radiation reaction 

In this category the radiation reaction is fully taken into account even after the contact. However 
the computation is time consuming because to estimate the radiation reaction potential two more 
Poisson e qua tions other than the Poisson equation for the Newtonian gravitational potential should 
be solved M || ||, ||, ^ |30j |lj || . 



In Fig. lj we show an typical example of such a simulation |2q| . In this case each mass of the 
binary is the same. We start the simulation from the equilibrium obtained in category A). Due to the 
loss of the angular momentum by gravitational waves the merging starts spontaneously. As merging 
proceeds, two arm spiral extends outward. The spiral arms become tightly bound later and becomes 
a disk. The final result is an almost axially symmetric central object and a disk around it. In Fig|l] 
the thick circle shows the Schwarzschild radius of the total mass. Although it is not possible to say 
the formation of a black hole in Newtonian simulations, we expect that the central object is a black 
hole. So the final destiny may be a black hole and a disk. This behavior is seen also in many other 
simulations in category B) and C). In Figj2|we show a typical wave pattern of the gravitational waves. 

In Fig.|, a simulation for different mass case is shown [p7f. In this case the lighter neutron star 
(left in Fig.3-a) is tidally disrupted by the heavier neutron star. We see only one arm . However the 
final object is more or less similar to the equal mass case, that is, a black hole and a disk. This final 
object may be relevant to the central engine of gamma ray burst. Ruffert will discuss in this volume 
on the coalescing binary neutron stars and gamma ray bursts using their simulations including the 
emission of neutrinos fl3C| [sT| . 

D) A)+1PN effect of Post Newtonian Hydrodynamics and Radiation Reaction 

Since neutron stars are general relativistic compact objects, the Newtonian hydrodynamics is 
zeroth approximation to the problem. One possible refinement is to include 1PN force [i.e. O ((u/c) 2 ) 
effect of the general relativistic gravity after the Newtonian gravity] . Oohara and Nakamura performed 
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Figure 1: Density and velocity on the x-y plane for EQSP2. The time in units of milliseconds is 
shown. Arrows indicate the velocity vectors of the matter. A fat line shows a circle of radius 2GM t /c 2 . 
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Figure 3: Density and velocity on the x-y plane for TD1. 



4 



Nowtonlan 



Post-Nowtonlan 



Nowtonlan Post-Nowtonian 




Figure 5: Wave forms of h+ and h x observed on the z-axis at lOMpc. The solid and dashed lines are 
for PN and N, respectively. 



such a simulation |32j. They first include 1PN force to the simulation in Fig.|l| and found that 
the 1PN effect is too strong. They decrease the mass of the neutron star to 0.62M Q froml.5M 
keeping the radius of the neutron the same so that the gravity becomes weaker. To see the difference 
between Newtonian and 1PN cases, in Fig]!] two simulations from the same initial data are shown 
J32| . We clearly found that they are different even in this rather weak gravity case. In Fig.| the 
wave forms are shown. The solid line (1PN) and the dashed line (Newton) are different. Due to the 
strong gravity effect of 1PN force, the strong shock in the central is formed and the coalescence is 
delayed. If the higher post Newtonian effects (2PN,...) are included, different results may be obtained. 
This suggests that fully general relativistic simulations are needed. This suggests also that all the 
results of coalescing binary neutron stars based on Newtonian hydrodynamics or Post-Newtonian 
hydrodynamics need to be viewed with caution. 

There is another point to be considered. Suppose the point particle limit of each neutron star. 
In the Newtonian gravity there exists a circular orbit for any small separation a. However in general 
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relativity a stable circular orbit does not exist for a < nsco wnere ISCO stands for Innermost Stable 
Circular Orbit. For a test particle motion in the Schwarzschild black hole nsco = QGM/c 2 where 
M is the mass of black hole. For equal mass point mass binary nsco is estimated as ~ 14Gm/c 2 
where m is each mass of the binary |Q. Since the radius of the neutron star is ~ 5Gm/c 2 , when 
two neutron stars contact, the separation a is ~ lOGm/c 2 , which is comparable to nsco m the point 
particle limit. Therefore in the last three milliseconds the general relativity is also important in the 
orbital motion. Recently ISCO for binary neutron stars including the effect of general relativity and 
the finite size of the neutron stars has been extensively studied. Shibata discusses ISCO problem in 
this volume. 



3 General Relativistic Simulations 

For full general relativistic 3D simulations of coalescing binary neutron stars, we need ISCO of binary 
neutron stars as initial data. However at present such initial data are not available. We therefore 
discuss the development of 3D general relativistic code and test simulations. We present only the 
summary of the basic equations. For details of formalisms, basic equations and numerical methods 
please refer |3~5fl . 

3.1 Basic Equations 
3.1.1 Initial Value Equations 

We adopt the (3+l)-formalism of the Einstein equation. Then the line element is written as 

ds 2 = -a 2 dt 2 + (dx l + f3 l dt) {dx 3 + 13° dt) , (5) 

where a, j3 l and 7y are the lapse function, the shift vector and the metric of 3-space, respectively. 
The lapse and shift represent the coordinate degree of freedom and its choice is essential for numerical 
simulations. 7^ is dynamical variable including gravitational waves. 

We assume the perfect fluid. The energy momentum tensor is given by 

T^v = {p + ps+p)u^u v +pg^ v , (6) 

where p, e and p are the proper mass density, the specific internal energy and the pressure, respectively, 
and Ufj, is the four-velocity of the fluid. We define p H = n tl n v T ilv , Ji = ~h i ^n v T illv and = 
hi^hj^T^ where n M is the unit normal to the i=constant hypersurface and = g^ v + n^n v is the 
projection tensor to the hypersurface. 

The initial data should satisfy the Hamiltonian and momentum constraint equations given by 

R + K 2 - Kij = 16tt Ph , (7) 
Dj (K j i - 8UK) = 8TT,J t , (8) 

where R is the 3-dimensional Ricci scalar curvature, is the extrinsic curvature and Di denotes 
the covariant derivative with respect to 7^ . 



3.1.2 Relativistic Hydrodynamics 

In order to obtain equations similar to the Newtonian hydrodynamics equations, we define p N = 
y/jau°p and — Ji/au°p where 7 = det(7^). Then the equation for the conservation of baryon 
number is expressed as 

d t p N +d e {p N V i )=0 1 (9) 
where V e — u /u°. The equation for momentum is 

dt{p N uf) + d e (p N ufV t ) = -^yad l p - ^(p + p H )dia 

^aJ k J e 



2(P + P H ) 



d llk i + • (10) 
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The equation for internal energy becomes 

d t (p N e) + d e (p N eV e ) = -pd v {^au") . (11) 
To complete hydrodynamics equations, we need an equation of state, p = p(e, p) 

3.1.3 Time Evolution of the Metric Tensor 

The evolution equation for the extrinsic curvature which is essentially the time derivative of 7y is 
given by 

d t Kij = a{R ij -8Tr[S ij + ^i j (p H -S e i )]}-D i D j a (12) 
+ a(KK ij -2K u K* j ), (13) 
+ K mi d 3 {5 m +K mj di0 m -K l3 d m (3 m , (14) 

where i£y is the 3-dimensional Ricci tensor. 

Now we define the conformal factor as (j) = (det [7ij]) 2 and 7^ as 7^ = 7^ /0 4 . The conformal 
factor <fi is determined by 

A0 = (l6np H + KijK ij -K 2 - , (15) 

where A and R is the Laplacian and the scalar curvature with respect to 7y, respectively. The time 
evolution of 7y is given by 

dfiij = -2a (itij - pijKj + Difa + DjJk - ^IijDiP 1 (16) 
= A T 

— ST-lj 1 

where 

K i6 =^K ih Pi = (fr^Pi (17) 
and Di denotes the covariant derivative with respect to 7^ . 

3.1.4 Coordinate Conditions 

On of the most serious technical problems in numerical relativity is the control of the space and 
time coordinates. Difficulties arise from (a) coordinate singularities caused by strong gravity and 
the dragging effect, where coordinate lines may be driven towards or away from each other, and (b) 
spacetime singularities which appear when a black hole is formed. We should demand coordinate 
conditions so that these singularities can be avoided. 

Spatial coordinates — the shift vector P l 



We demand diAf } — 0. Then Af } becomes transverse-traceless. Equation (|l(j) is reduced to the 
equation for the shift vector P % as 

v 2 p l + ±di (d e p e ) = 



9 j 



2a (ku - p i3 K 

2 

%$iP + %djP e - -%d £ p e + P i dt% 



(18) 



where V 2 is the simple fiat-space Laplacian and 7^ = 7y — 5 



13 ■ 
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Time Slicing — the lapse function a 



We require a time slice to have the following properties: (1) In the central region, spacetime 
singularities should be avoided. These singularities are regions of spacetime where curvature, density 
and other quantities are infinite. (2) In the exterior vacuum region, the metric should be stationary 
except for the wave parts to catch gravitational waves numerically. 

The lapse function a is given by 



a = exp 




(19) 



where (j) — (j) — 1. In this slicing, the space outside the central matter quickly approaches the 
Schwarzschild metric. Since the non-wave parts decrease as 0(r -4 ) for large r, we can catch grav- 
itational waves numerically with 0((M/r c ) 3 ) accuracy, where M is mass of the system and r c is 
distance from the origin of the position where gravitational waves are evaluated. 



4 Coalescing Binary Neutron Stars in 3D Numerical Relativ- 
ity 

We show our test simulations of coalescing binary neutron stars. We place two spherical neutron 
stars of mass M = I.OMq and radius ro = QM with density distribution (p N ) of 7 = 2 polytrope at 
x = 7.2M and x — —7.2M. We add a rigid rotation with angular velocity SI as well as an approaching 
velocity to this system such that 

f —ty - if x > 0, 

U * = o n •< n (20) 

[ -Sly + Slr if x < 0, 

v£ = fix. (21) 

We adopt (x,y,z) coordinates with 201x201x201 grid size (9GB) and typical CPU time being 10 
hours on VPP300(FUJITU). The details of a numerical method is shown in ]35| . 

We performed two simulations, BI1 and BI2, with different values of SI. The total angular momen- 
tum divided by the square of the gravitational mass are 0.99 and 1.4 and for BI1 and BI2, respectively. 
Figure 6 shows the evolution of the density on the x-y, y-z and x-z planes of BI2. The binary starts to 
coalesce and rotates about 180 degree. In Fig.^ we show flux = r 2 ^ i j/32n which can be considered 
as the energy flux of the gravitational wave energy if we take appropriate averaging, flux shows an 
interesting feature of generation and propagation of the waves. A spiral pattern appears on the x-y 
plane while different patterns with peaks around z-axis appear on the x-z and y-z planes. This can 
be explained by the quadrupole wave pattern. In Fig.|^ we show the luminosity of the gravitational 
waves as a function of the retarded time t — r calculated at various r. For large r it seems that 
the luminosity converges, which suggests that the gravitational waves are calculated correctly in our 
simulations. 

The total energy of the gravitational radiation as a function of the retarded time is shown in Fig.^J. 
The total energy amounts to 4 x 10~ 3 M (0.2% of the total mass) and 1 x 10 _3 M (0.05% of the total 
mass) for BI2 and BI1, respectively. 

In conclusion, the present test simulations suggest future fruitful development of 3D fully general 
relativistic simulations of coalescing binary neutron stars. 

Numerical simulations were performed by VPP300 at NAO. This work was also supported by a 
Grant-in- Aid for Basic Research of the Ministry of Education, Culture, and Sports N0.O8NPO8OI. 
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Figure 6: Mass density pjy on the x — y 
plane. Arrows indicate the velocity vector 
V . Quantities are represented in units of 
c = G = M Q = 1. 
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Figure 7: r 2 7 i:) /327r, which is expected to rep- 
resent "Energy density of the gravitational 
waves." 
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Figure 8: Luminosity of gravitational waves Lgw 
as a function of the retarded time t — r calculated 
at r = 30 - 50. 
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Fi gure 9: Energy of gravitational waves i?GW ^ 
a function of the retarded time t — r calculated 
at r = 30 - 50. 
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